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Computer simulations in QCD are based on the discretization of the theory on a Euclidean lattice. 
To compute the mean value of an observable, usually the Hybrid Monte Carlo method is applied. 
Here equations of motion, derived from an Hamiltonian, have to be solved numerically. Com- 
monly the Leapfrog (Stoermer-Verlet) method or splitting methods with multiple timescales a la 
Sexton- Weingarten are used to integrate the dynamical system, defined on a Lie group. 
Here we formulate time-reversible higher order integrators based on implicit partitioned Runge- 
t— \ Kutta schemes and show that they allow for larger step-sizes than the Leapfrog method. Since 

these methods are based on an infinite series of exponential functions, we concentrate on the 
truncation of this series with respect to the global error and accuracy. Finally, we see that the 
global error of a SPRK scheme is always even such that a convergence order of one is gained for 
methods with odd convergence order. 
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1. Introduction and Motivation 

In the molecular dynamics step of the Hybrid Monte Carlo method [1], Hamiltonian equa- 
tions of motions have to be solved. These equations form coupled systems of matrix differential 
equations of the form 

0r= 3HMM (Ua) 

oA v 

dH([U], [A]) , , 
A v = at/ =g(^v). farv = l,...,n. (1.1b) 

In this notation, U v is an element of a matrix Lie group G and A v an element of its associated Lie 
algebra q. Thus, [U ] can be imagined as a vector of n matrix Lie group elements U\,U2,..., U n , and 
[A] as a vector of n Lie algebra elements A\,A2, ■ ■ ■ ,A n . In pure lattice gauge theory, the element 
U v can be seen as the link U xfl between the lattice sites x and x + ajX. Thus, A^ is its associated 
momentum P x ^ times the complex i. In this context, the vectors of matrices [U ] and [A] are the 
whole configurations of the links and its momenta. 

The equations of motion have to be solved in a Lie group, respectively in a Lie algebra with a time- 
reversible and area-preserving scheme. In a recent paper [2], we have investigated the potentialities 
of higher order partitioned Runge-Kutta schemes for solving the equations of motions such that the 
desired properties are met. We found out that symmetric partitioned Runge-Kutta methods based 
on the Magnus and Munthe-Kaas approach can be time-reversible. So far, area-preservation is not 
fulfilled and must be corrected in the acceptance step. Furthermore, the global error of this scheme 
is always even and investigated in detail in this paper. In doing so, we start with a short derivation 
of symmetric partitioned Runge-Kutta schemes based on the ideas of Magnus and Munthe-Kaas. 
Afterwards, we focus on the global error and accuracy of the method and show some numerical 
results. 



2. Numerical Integration 

The differential equations (1.1) become an initial value problem (IVP) by prescribed initial 
values: U v (0) := £/ v .o an d A v (0) := A v .o for v = 1, . . . ,n. Thereby, the initial values U v (0) have 
to be in the Lie group and the elements A v (0) in the Lie algebra. Considering the structure of 
equation (1.1), (1.1a) is a differential equation in a Lie group such that it has to be solved with a 
numerical scheme that guarantees a solution in the Lie group as described in paragraph 2. 1 For the 
second equation (Lib), no special treatment has to be applied. It is an equation in the Lie algebra 
0, which is a linear space. Thus, this equation can be solved with any time-reversible and area- 
preserving numerical scheme. For convenience, we leave out the index v from now on. This means 
we investigate just one coupled differential equation for a special but arbitrary index v: U =A-U 
and A = g(U). The results can then be extended straightforward to the whole vectors [U] and [A]. 

2.1 Differential Equations in Lie Groups 

Concerning equation (1.1a), we follow the ideas of Magnus and Munthe-Kaas. Magnus [3] 
stated that the differential equation (1.1a) in the Lie group can be replaced through a differential 
equation in the Lie algebra. This new differential equation can be solved directly due to the linearity 
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of the Lie algebra. The strategy is as follows: Identify U (t) with exp(0(/)) such that the variable 
changes from U to £2. £2 is the solution of the differential equation 

D. = dexp n l (A), (2.1) 

with Cl(t) G fj and initial value £2(0) := 0. The derivative of the inverse exponential map (2.1) is 
given by an infinite series as 

Jexp n 1 (A) = £^a4(A). 
k>0 K - 

In this series, the variables are the k-th Bernoulli numbers and the adjoint operator ad^ is a 
mapping in the Lie algebra q given by ad^(A) := [£2,A] = £2A — A£2. It follows the conventions 
ad£(A) =A and ad£>(A) = [fl,ad^ _1 (A)]. This means, £2 is the solution of the differential equation 

k>0 K - 

Knowing £2, the solution U of (1.1a) can be attained via U = exp(£2)£/o- In total, we record that 
the initial value problem (1.1) is equivalent to 

£2 = £ ad^(A), A = g(U) with [/ = exp(£2)£/ (2-2) 

/t=o 

and f/ (0) := t/o G G, A(0) := Ao G and £2(0) := G 0. This transformed problem can now be 
solved directly by a Runge-Kutta method without destroying the Lie group structure: As the Lie 
algebra is a vector space, the analytic solution (£2(7), A(f)) as well as its approximation (£2i,Ai) 
attained by a numerical integration scheme both are elements of the Lie algebra 0. Furthermore, as 
for any a G the matrix exponential exp(a) is in the associated matrix Lie group G, also U is in G. 

2.2 Symmetric Partitioned Runge-Kutta schemes 

The problem in solving (2.2) is that £2 is given as infinite series which has to be suitably 
truncated after q + 1 terms. This means, the truncation index q of £2 has to be chosen properly 
such that a numerical integration scheme meets a prescribed convergence order p. Thereby, the 
convergence order of a numerical integration method is p if the deviation between the exact solution 
and its numerical approximation after one step is of order p + 1 in a suitable norm. Here, the idea 
of Munthe-Kaas comes into play. He states in [4] that the truncation index q of £2 has to be chosen 
as a value larger than the desired convergence order p minus one. Consequently, for a Runge-Kutta 
scheme of convergence order p, £2 is set as a function depending on the truncation q > p — 2 of the 
aforementioned infinite series, i. e. 

h = i,Tr ad $ x (A)=:f q (Q,A). (2.3) 
k=o K - 

All in all, the exact solution of (2.2) is approximated through an integration scheme of order p of 
the truncated model 

« = *E" ff^fi^)' A=g(U) with U = exp(£2)£/ , (2.4) 
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0(0) := Uo G G, A(0) := A G and H(0) := G 0. Thereby, this model can be solved with 
higher order time-reversible symmetric partitioned Runge-Kutta (SPRK) schemes derived in [2] as 
follows: Compute the approximations 

a^h^biKi, A^Ao + h^biLt, (2.5) 

i=l i=l 

with increments Kj = f q (Q.j,Ai) and L, = g (Z?,-) for i = 1, . . . ,s. In the course of this, the internal 
stages are defined as 

s s s 

D. i = hJ^a ij K j , A i =A + hY,oc ij L j , U i = exp(X i )exp(±ei l )U , X i = hY d YijKj- 

;'=i j=i j=i 

At the end, the solution U\ is attained via U\ = exp(Q.i)Uo. In this scheme, the coefficients 
bi,bi, CCij, CCij and Yij f° r hj = 1, • • • , J can be determined to guarantee time-reversibility (and sym- 
metry). Their values for convergence order p = 3 can be found in [2]. 

3. Global Error and Accuracy of the SPRK Method 

For the local error, the solution of the integration method after one step has to be compared 
with the exact solution U(to + h),A(to +h) of the differential equations (1.1). The SPRK method 
(2.5) is of convergence order p if 

\\U(t + h)-Ui\\ = 0(h p+i ) and \\A(t + h) -A t \\ = 0(h p+l ) (3.1) 

holds. Since the approximation to U\ is computed from evaluating the matrix exponential (we 
assume here that we can evaluate this exactly), which is Lipschitz on every closed interval, it 
suffices to demand 

\\a(tQ + h)-Q. x \\ = G(h p+V ) and ||A(f + h) -A x \\ = G(h p+X ) 

with exact solution £l(to + h),A(to + h). According to Munthe-Kaas, the approximations £2i and 
Ai of the exact solution of the suitably truncated problem (2.4) can also be interpreted as approxi- 
mations to the original problem (2.2). With the same argument, we can even formulate a stronger 
statement on the local accuracy (3.1). As the method is symmetric, theorem 3.2 in [6, II.3] applies, 
which states that the maximal convergence order p of a symmetric method is even, which means 
that the local error is always odd. Hence, the SPRK method developed as a method of an odd con- 
vergence order p is of order p + 1 . The global error of a numerical integrated scheme is computed 
as the sum of the local errors. This means, for the computation of a trajectory with length z, an 
integration method with fixed step size h is applied N = zjh times. Hence, the global error is of the 
order local error minus one. Thus, the SPRK method of convergence order p has at least a global 
error of order p. Again, in case of an odd p, the global error is of order p + l. 
The accuracy depends on the truncation k = q = p — 2 of the series given in (2.2) with p being 
the convergence order of the Runge-Kutta method that is used to solve the problem numerically. 
We recall the basic steps of the proof given in [6] for a deeper understanding of the choice of the 
truncation parameter q. For this purpose, we restrict to an uncoupled Lie algebra problem 

n = F 00 (Q):=£-^ad^(A), with £1(0) = G 0, (3.2) 
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which arises from the Magnus approach mentioned in paragraph 2. 1 . The truncation of the series 
in (3.2) at k = q yields the truncated problem 



such that 



Q = F q (Q), with £1(0) = 0, 



F.(fl(0)-F f (n(0)= I ^< (r) (A). 



(3.3) 



k=q+l 



For sufficiently smooth A we recognize 



Sl(t) - Q(t) || = ||fl(0)+ / F„(Q(T))dT-(fl(0)+ / F q (a(r))dr 



F 00 (0(T))dT- 



F 00 (n(T)) + C(T)dT 



ad^ ) (A) = ^+ 1 ), 
which is due to the nested structure of the ad-Operator [6]. Hence, we have 

F q (h(t))=F BO (Q(t))+C(t) 

with C{t) = c\t q+2 + C2? ?+3 H and constant values c\,c<i, ■ ■ ■ G 0. For fixed h > and f G [0,h] 

the exact solutions of (3.2) and Q.(t) of (3.3) satisfy 



(3.4) 



The function F M is Lipschitz continuous on every closed interval for sufficiently smooth A. We 
assume that for an interval where both Cl(t) and £2(?) reside in for t G [0,h], the Lipschitz constant 
is Loo G M, i. e., 



||F oo (n(0)-F»(n(0)|| < M n (0 -^(Oil- 
Furthermore, for ? G [0,/i] we see that 

||C(0ll < ||ci||/j 9+2 + ||c 2 ||/i' ?+3 + ... :=c(h)EM. 
Hence from (3.4) it follows that 

||n(*)-n(*)|| <c(h)-t+Loo [ ||n(T)-S(T)||dT, 



< / ||F»(fl(T)) -F»(Q(t)) II dT + / ||C(T)||dT 



(3.5) 



such that the requirements of the "Gronwall lemma" [5] are satisfied by which 

m) -a(t)\\<^(e L °° h -i 



±- (||d ||/^+ 2 + ||c 2 ||//' +3 + •••)•((! + ^ + ^(Z^/i) 2 - 
(||d w +1 + ||c 2 ||^+ 3 + •••)• (h + ' L M /* 2 + ^L^ 3 + • • 
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Thus, the difference between the exact solutions of the full problem (3.2) and the truncated problem 
(3.3) is 

\\£l(h)-&{h)\\ = &(h q+i ) (3.6) 

after one time step h. Applying a one step method of convergence order p on the truncated problem 
(3.3) means to calculate an approximation Q.\ to the exact value £l(h) such that 

\\h(h)-ai\\ = @(h p+l ). (3.7) 

Finally, we interpret Q.i as an approximation to the exact solution of the original problem (3.2). 
The quality of this approximation is determined by the deviation (3.6) introduced by the modeling 
and the discretization error (3.7): 

||n(A)-£2i|| < ||n(fc)-Q(fc)|| + ||Q(ft)-fli|| = @{h q+3 ) + @{h p+l ). 

This clearly indicates that Q.i is a numerical approximation to Q.(h) of convergence order p, i. e., 

= <ff(h p+l ) if g + 3 >/> + !, i.e.,q>p-2. 



4. Numerical Tests 

We consider a pure lattice gauge theory in SU(2,C) with Wilson action and compare the SPRK 
method described in (2.5) with the Leapfrog method. For this purpose, we investigate a symmetric 
partitioned Runge-Kutta scheme of convergence order p = 4 which contains the truncated function 
Cl = fq(Q.,A) given in (2.3). Because of the symmetry, the method has an even convergence order 
such that the choice p = 3,q = 1 already leads to a local error of order 5. This means, we use the 
equation 

d = f 1 (a,A)=A-^[a,A] (4.i) 

according to (2.3) and perform simulations on a 2-dimensional lattice with lattice size L = T = 32. 
There are 2 results shown in figure 1: On the left side, the convergence order of the different 
methods can be seen. For this purpose, we consider the energy change AH of two successive con- 
figurations after a whole trajectory of length 1 and take the mean of 5000 configurations. The 




Figure 1: Left: Convergence order. Right: Area-preservation. 
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statistical errors are so small that they are not visible in the plot. Since the energy change AH de- 
viates from zero just because of the numerical errors of the integration scheme, the violation of the 
energy preservation gives the global error. We see that for a given \AH\ the SPRK allows for larger 
step sizes. On the right side of figure 1, we see the violation of the area-preservation in dependence 
of the step sizes chosen in the numerical methods. Area-preservation (up to roundoff errors) is 
given if the determinant of d(D.i,A\)/d(D^,Ao) has exactly the value 1. Here, the determinant is 
numerically approximated by first order difference quotients. 

5. Conclusion 

We investigated the accuracy of the time-reversible symmetric partitioned Runge-Kutta scheme 
(2.5). The order of accuracy consists of two components: On the one hand, the convergence order 
depends of course on the order p of the method itself. As the method is symmetric, the local error 
is always odd, i. e., the scheme has a local error of order p + 1 for an even convergence order p. On 
the other hand, the SPRK scheme contains one truncated series (2.3). The truncation index q has 
to be larger than or equal to p — 2 to meet the prescribed convergence order. All in all, choosing 
an SPRK method with an even local error should be preferred since a convergence order of one is 
gained by the symmetry. We performed simulations for an SPRK scheme of convergence order 4 
and see that the global error given in the numerical results has order 4 as theoretically expected. 
In the development of the SPRK method, area-preservation has not been considered. Thus, it is 
not surprising, that area-preservation is not met applying this scheme. This property has to be 
investigated in future work. 
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